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The role of inhibitory feedback for information processing in thalamocortical circuits 
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The information transfer in the thalamus is blocked dynamically during sleep, in conjunction with 
the occurrence of spindle waves. In order to describe the dynamic mechanisms which control the 
sensory transfer of information, it is necessary to have a qualitative model for the response properties 
of thalamic neurons. As the theoretical understanding of the mechanism remains incomplete, we 
analyze two modeling approaches for a recent experiment by Le Masson et al. on the thalamocortical 
loop. We use a conductance based model in order to motivate an extension of the Hindmarsh-Rose 
model, which mimics experimental observations of Le Masson et al. Typically, thalamic neurons 
posses two different firing modes, depending on their membrane potential. At depolarized potentials, 
the cells fire in a single spike mode and relay synaptic inputs in a one-to-one manner to the cortex. 
If the cell gets hyperpolarized, T-type calcium currents generate burst-mode firing which leads to 
a decrease in the spike transfer. In thalamocortical circuits, the cell-membrane gets hyperpolarized 
by recurrent inhibitory feedback loops. In the case of reciprocally coupled excitatory and inhibitory 
neurons, inhibitory feedback leads to metastable self-sustained oscillations, which mask the incoming 
input, and thereby reduce the information transfer significantly. 

PACS numbers: 87.19.La 05.45.-a 87.19.Nn 84.35.+i 



I. INTRODUCTION 

Spindle Oscillations are waxing and waning waves that 
originate in the thalamus at the transition to slow wave 
sleep in mammals Q. In a pioneering experiment 0, 
spindle oscillations have been investigated in a hybrid 
thalamic circuit consisting of a biological thalamocorti- 
cal (TC) and an artificial reticular (RE) thalamic cell . 
An artificial cell was necessary to firstly manipulate selec- 
tively the synaptic connection of the inhibitory feedback 
from the RE to the TC cell (Fig-QJ. Spindle oscillations 
depended critically on the presence of the synaptic con- 
nections between TC and RE cells 0,0. Further, the 
occurrence of spindle oscillations comes along with a sig- 
nificant decorrelation between input and the output of 
the TC cells. Clarification of the relationship between 
the occurrence of spindle oscillations and the decrease in 
information transmission could help to gain more insight 
into the mechanisms which deprive the sensory informa- 
tion from the consciousness while mammals are sleeping. 

Spindle, or bursting, oscillations in excitable media 
are widespreadly observed in a variety of physical [5j, 
chemical and biological oscillators. Typically, the 
potential of a bursting neural cell undergoes subsequent 
shifting between active and silent phases. In the active 
phase, the membrane potential oscillates quickly, and in 
the silent phase, it evolves slowly without or only with 
subthreshold oscillations 9j. In this work, we focus on 
bursting or spindle oscillations in neural systems, and 
investigate their role for information transfer through 
the thalamus. The thalamus, our major gateway for all 
sensory information, relays the incoming information de- 
pending on the state of arousal jToL fill . To process and 
relay information |12| is a generally important, but not 
equally well understood, feature of neural systems; here 
we contribute to the understanding of the underlying dy- 



namics of a key mechanism in the thalamocortical loop. 
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FIG. 1: Structure of synaptic connections in a pair of recip- 
rocally coupled RE-TC cells. 



The paper is organized as follows. In section [HI we 
use biophysical models of TC and RE neurons intercon- 
nected with realistic model synapses as introduced in Q 
and also used in the hybrid network investigated by Le 
Masson in 0- In sect ion ILTT] we establish an extension of 
the Hindmarsh-Rose model jl4| . For this purpose, we use 
dynamically coupled Hindmarsh-Rose neurons to repro- 
duce the experimental results of 0- We show that the 
Hindmarsh-Rose system has to be extended by a fourth 
variable working on a very slow timescale, in order to 
describe the experimental results. 

For both models, we investigate the influence of the 
inhibitory synaptic connection on the information trans- 
fer of the TC cell and compare the results with the ex- 
periment of Le Masson et al. 0- An additional goal 
is to answer the question whether the complicated bio- 
physical conductance-based model can be replaced by the 
much simpler Hindmarsh-Rose model without loss of the 
most important dynamical features. We show how to 
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achieve this by introducing an extended Hindmarsh-Rose 
model with an additional degree of freedom. Some im- 
portant features are better reproduced by the extended 
Hindmarsh-Rose neuron than by the biophysical model. 
Our model is by far less computationally costly, yet sim- 
ple enough to understand dynamical mechanisms. 



II. INFORMATION TRANSFER IN A 

BIOPHYSICAL MODEL OF A 
THALAMOCORTICAL OSCILLATOR 

A common approach to model the thalamocortical sys- 
tem are conductance-based models EE EH)]. De- 
scriptions of the properties of single neural cells can be 
found for example in |3|, |J| . The case of two intercon- 
nected cells without external forcing has been investi- 
gated in £| . To our knowledge the case of two intercon- 
nected cells as shown in Fig. ^ which get excited by an 
unregularly distributed realistic synaptic bombardment 
as used in the experiment of Le Masson et al. has not 
been studied in detail up to now. In this section we will 
investigate the influence of the inhibitory synaptic con- 
nection on the information transfer of the TC cell for this 
case. Destexhe et al. model the TC cell and the RE 
thalamic cell by using a conductance-based single com- 
partment model. The time evolution of the membrane 
potentials is governed by the following cable equations 

C m Vr = —It l — It — hi — It Ncl — It k 

-lGABA a - IcABA b (1) 

C m Vn = —Ir l — Its — Ik[Co\ — Ican — In. Na — Ir k 
-Iglu, (2) 

where Vt is the membrane potential of the TC cell; C m 
is the capacity of the membrane. According to 0, we 
included the following ion currents: The leakage cur- 
rent It l and Ir l , the low-threshold Ca 2+ currents It 
and Its, the hyperpolarization activated current Ih, the 
Ca 2+ -activated currents Ik[Co\ and Ican, and, like in 
the Traub and Miles E model, the fast Na + and K+ 
currents It/b. Ncl and It/r k , which are responsible for the 
generation of action potentials. The synaptic currents 
IcABA a and lGABA b represent the GABA a and GABAb 
receptors in the synapses from RE to TC cells, while 
Iglu describes the excitatory synapse from the TC to 
RE cells. 

Destexhe et al. |4| refer to several sources for the de- 
scription of the ion currents. The kinetic equations of 
the TC cell have been described in detail in El- The 
It kinetics was taken from Wang et al. l35|. Ih was 
described by the model of Destexhe et al. El> which 
incorporates a regulation by intracellular calcium. For 
the RE cell Its was taken from Huguenard and Prince 
El- The kinetics of the Ca 2+ -dependent currents Ix[Ca] 
and Ican are adjusted to the clamp data of RE neurons 
[l9| . For both cells, I^a and Ik were taken from |28|. 
The calcium dynamics in both cells are described by a 



simple model which was introduced in [lg. All details 
are described in appendix 1X1 

In biology, the typical input signal for the TC cells 
is a spike train with unregularly distributed inter-spike 
intervals, which can be modeled by a Poisson process 
with a refractory period. For example, the retinal cells 
fire such unregularly distributed spikes in darkness or 
under constant illumination. The refractory period is 
necessary because every biological cell needs some time 
r r to recover after it has fired a spike, so the cell cannot 
fire spikes with arbitrary low inter spike intervals. In our 
model, the ISI are distributed according to the following 
distribution 







for r < T r 

for T > Tr 



(3) 



where e rrr is a scaling factor, f = 1/100 and r r = 30ms 
here. If we stimulate our model by such a modified Pois- 
son process, we have a computational model for the ex- 
periment of Le Masson et al. which allows us to 
compare our computational model with the experiment. 



A. Simulation results of the biophysical model: 
Spindle oscillations, waxing and waning, and 
influence on information transfer 

With the TC and RE cells interconnected, the compu- 
tational thalamic circuit receiving this realistic synaptic 
bombardment Q showed waxing and waning oscillations. 
These spindle sequences consisted of bursts of 8-10Hz os- 
cillations lasting a few seconds, and were separated by 
silent phases of around 8-10s. The periodical occurrence 
of this oscillations is very similar to biological spindle 
waves which occur in sleep-related states. As in the ex- 
periment, each spindling state is terminated by an after- 
depolarization of the membrane potential. The reason 
for this effect is the hyperpolarization activated current 
Ih (see appendix Ia"|| ; when this after-depolarization adds 
to the synaptic bombardment, it forms a depolarization 
which deactivates the low threshold Ca 2+ current It and 
subsequently blocks the input signal transfer. 



B. Measuring the information transfer 

During the spindling state, the firing pattern of the TC 
cells (excitor) , which is very different from the input (see 
Fig. Et~b), shows that obviously the information trans- 
fer of the input is low. Le Masson et al. Q answer the 
question whether this low transfer is still reliable, by cal- 
culating two different indices. First, the contribution in- 
dex (Tci) examines the TC cell discharge and quantifies 
the percentage of output spikes N out which are precisely 
correlated with retinal input spikes. It estimates the reli- 
ability with which a TC spike can be considered as being 
triggered by an input spike rather than being sponta- 
neous. It is computed as the peak of the crosscorrelation, 
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FIG. 2: (Color online) (a) Spontaneous spindle activity in 
the computational thalamic circuit. The TC cell receives ar- 
tificial synaptic retinal bombardment, modeled by a Poisson- 
distributed spike train. The RE-TC oscillator shows spindle 
activity, (b) Detail of (a) . (c) Total synaptic current injected 
into the biological TC cell and simultaneous voltage of the 
TC and RE cell. Bottom: The hyperpolarization of the TC 
cell activates the It current, what initiates a postinhibitory 
rebound burst (see also Fig.|SJ. Calibration bars: 10s, 20mV 
(a); Is, 20mV (b); 10ms, 6nA, 20mV, InA (c). 



normalized by the number of output spikes. Second, the 
correlation index {Tec) measures the global efficiency of 
the input-output spike transfer and indicates the ratio of 
input spikes being actually transmitted as output spikes 
in the TC neuron. It is computed as the peak of the 
crosscorrclation between retinal and TC neuron spikes, 
normalized by the number of retinal cell spikes. To be 
able to compare our numerical results with the experi- 
mental results of Le Masson et al. Q quantitatively, we 
will use exactly these measures here. 

To be able to compare the numerical results obtained 
from the biophysical model with the simplified model in- 
troduced in this paper, we will define two new measures 
for the information transfer which are equivalent to the 
measures used by Le Masson et al. in 2], but do not 
depend on the detailed form of the signals. First, the 
signal-to-noise ratio (Tsn) examines the amplitude of the 
TC-cell signals and quantifies the percentage N tr of out- 
put spikes N out which are exactly triggered by an input 
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FIG. 3: (Color online) Averaged contribution {Tci) (a) and 
correlation {Tec) (b) indices versus synaptic strength. 



spike; such it is a measure equivalent to the Tci used. 
An output spike is considered to be triggered by an in- 
put spike if it occurs within a delay of < 50ms after an 
input spike. As the Tec-, it is a measure for the reliability 
with which an output spike can be considered as being 
triggered by an input spike rather than being due to the 
intrinsic dynamic of the circuit. 



Tsn — 



N,, 



(4) 



Second, the transfer efficiency Tte measures the global 
efficiency in the input-output spike transfer, and indi- 
cates the ratio of input spikes being actually transmitted 
as output spikes, obviously it is equivalent to the measure 
Tec- It is computed by the number of triggered spikes 
N tr divided by the total number of input spikes, 



Tte — 



Nt, 
N„ 



(5) 



The difference between the measures (Tsn) and {Tci) is 
that our measure checks if an output spike occurs within 
a certain delay after an input spike, while the contribu- 
tion index counts the spikes at a fixed time delay defined 
by the maximum of the crosscorrelation. The Tec and 
Tte differ in the same way as mentioned above, further 
Tte only takes values between zero and one. As our 
measures only count the spikes, and such treat spikes like 
binary events, they can be used for any spiking system 
and allow a comparison of the transmission properties of 
different models. The disadvantage is that they do not 
shed light on the amplified details of a signal system such 
as Tci and Tec do. 

These two measures, the transfer efficiency and the 
transfer reliability, are necessary to characterize the in- 
formation transfer appropriately. The usage of two dif- 
ferent transfer measures allows to answer the question 
whether the transfer gets blocked or if incoming signals 
get masqued by autonomous oscillations of the system. 
The results are given in Fig. 
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In the presence of strong inhibitory feedback 
{gaABAa = O.lmS and g GA BA b = O.OlmS) the (T CI ) and 
(Tsn) were low (« 0.4), indicating that less than half of 
the output spikes were triggered by an input spike; thus 
the thalamus is not transferring spikes in a one-to-one 
manner. To answer the question whether the degree of in- 
hibition produced by the RE neurons could directly con- 
trol the precision of spike transfer, we varied the strength 
of inhibition. Our simulations show that the influence of 
the inhibitory feedback is by far not so strong as in the 
experiment [2J. In our computational model, the Tsn 
varied only between rs 0.6 and ss 0.4. In agreement with 
the experimental data 0, the global efficiency of spike 
transfer Tt e was not significantly different in the absence 
or presence of inhibitory feedback. 
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FIG. 4: (Color online) Percentage of output spikes triggered 
by an input spike (Tsn) (solid) and Transfer Efficiency Tte 
(dashed) for different values of the inhibitory feedback. Both 
indices show no significant change when the inhibitory feed- 
back is varied within the stability borders. 



C. Discussion of the biophysical model 

Altogether, our simulation cannot reproduce the trans- 
mission behavior of the thalamic circuit observed Q ex- 
perimentally. When we compare the numerical achieved 
values for the Tqi and Tec (Fig- El quantitatively with 
the values from the experiment, we find that the effect 
of input-output spike decorrelation is present, but lacks 
some features. First, even without or very low inhibitory 
feedback, the Tec is not more than w 0.25; in the ex- 
periment the Tec varied in a range from about 85-45 
(see Figure 3e in [2] ) . For a biological system this would 
mean that the animal would always be drowsy. Second, 
the system does not allow precise up and down regu- 
lation of the signal reliability in a wide range, the Tci 
only varies in range from around 0.7-0.35, that is about 
a factor 2, while in the experiment the Tqi varies from 
about 0.8-0.2 (see Figure 3d in 0) which corresponds 



to a factor 4. Further, spindle oscillations occurred even 
with very little or no inhibitory feedback, which is also 
in contrast to the experimental results. 

At this point it appears questionable whether the de- 
tailed biophysical model described above — with its high 
dimensionality and parameter space quite too detailed 
to gain systematic understanding — provides the appro- 
priate level of description to analyze coupled systems, 
e.g. of the thalamocortical loop investigated here. We 
will however use the biophysical model to establish our 
simplified model. As the biophysical model incorporates 
many features of realistic biological neurons, it allows us 
to identify the key elements which may be responsible 
for sleep spindle oscillations. The slow repetition rate of 
spindling (0.3-0.1 Hz) is due to intrinsic mechanisms of 
the TC cell. The s low va riables F 2 and S 2 (Fig. [TTJ) of 
the Ih current (see IIII Dl and appendix A) play the key 
role for these long silent phases. The spindling starts in 
the TC cell due to the hyperpolarization low threshold 
Ca 2+ current It (see Fig. and Fig.|3J. The spindling 
terminates due to the activation of the Ca 2+ activated 
current Ih (see Fig. EJ. The frequency of the spindling- 
is determined by intrinsic properties of the RE cell, as 
shown by Destexhe in but also the intrinsic proper- 
ties of the inhibitory synapses play a role. To provide 
an overview, Fig. [S] shows a schematic diagram of the 
mechanisms leading to spindle oscillations. 
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FIG. 5: Schematic diagram illustrating the generation mecha- 
nism of spindle oscillations. We distinguish two cycles: First, 
the spindle cycle is responsible for the spiking; the system runs 
several time through it before the waxing and waning cycle 
is activated, leading to bursts of oscillations. The slowly de- 
caying current Ih leads to the long silent epochs during which 
spike transfer is blocked. Altogether, the interplay of the cur- 
rents It and Ih plays an important role for the genesis of 
spindle oscillations, (see also Figs. and 0J. 



Unfortunately, due to the complexity of the high- 
dimensional biophysical system, the connection between 
biophysics and dynamics is hard to unravel. Further, the 
computational effort is very high, making studies of net- 
works almost impossible. To solve these problems, we 
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propose a reduced model to reproduce the experimental 
data of Le Masson et al. by using the simpler Hindmarsh- 
Rose type neurons, and modify them in a way that 
the features described above are reproduced. 

III. THE REDUCED SYSTEM: AN EXTENDED 
HINDMARSH-ROSE MODEL 

A. The Hindmarsh-Rose model 

The Hindmarsh-Rose equations are a simple polyno- 
mial model of bursting in thalamic cells, which reproduce 
several features, like for example rebound bursts, of more 
complicated biophysical models 0, 0, 13 . The original 
Hindmarsh-Rose equations are given by 

v = w - v 3 + 3v 2 - z + I(t) (6) 
ib = 1.8 — 5v 2 — w (7) 
z = e(3.3(u + 1.56) - z), (8) 

where v is the voltage, w is a recovery variable, I(t) is 
an external forcing and z is a slow variable. In order to 
understand the behavior of equations > we nrs t con- 
sider the limit e — > 0, so z(t) = z, which results in a 
planar reaction-diffusion system. This fast subsystem is 
responsible for the spiking dynamics of the Hindmarsh- 
Rose system. The dynamics of this subsystem is also 
crucial for e ^ 0. This procedure is called singular ap- 
proximation, or slaving principle, and is widely used to 
separate fast and slow subsystems [22, H3, H3 ■ Slow vari- 
ables (here z) can be treated as slowly varying parameters 
and the rest of the system can be studied as a function 
of these new parameters. We apply this method to the 
Hindmarsh-Rose model. 

To understand the influence of the variable z on the 
(v. w) subsystem, equations EHJl with I(t) = are trans- 
formed to a Lienard system piyj' yielding (Add. lB~2f) 

= v + f(v)v + g(v), (9) 

where g(v) — v 3 + 2v 2 — 1 can be considered as the gra- 
dient of a potential <!>(«, z), and f(v) = 1 — 6v + 3w 2 as a 
damping term. The roots of g(v , z) give the fixed points 
of the equation. Depending on the parameter z, there 
exist either one or three real roots. 

The bifurcation diagram of equations (JSJ-EI) as a func- 
tion of z is shown in Fig. |B| (the bifurcation diagram was 
computed with the auto interface [25| of xpp 26] ) . If we 
start with z = 2 from the fixed point (v = —2.05, w = 
— 19.67) of the fast subsystem and go to the left in the 
bifurcation diagram, the stable fixed point loses its sta- 
bility by a saddle node bifurcation at z w 0.6148. For 
values of z between [0.6148,0.8891], a stable limit cycle 
and a stable fixed point coexist. At z = —10.4, the up- 
per branch of the steady state loses stability by a Hopf 
bifurcation (for further details see 0,E3)- Here we will 
only focus on values of z between [—5, +5]. 
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FIG. 6: Maximum/Minimum bifurcation diagram of v in the 
fast subsystem (4) as a function of z. Heavy lines indicate sta- 
ble fixed points (limit cycles) and thin lines indicate unstable 
fixed points (limit cycles). (See Appendix IB 31 ) 



B. Periodic forcing of the Hindmarsh-Rose neuron 

Another important case is that of z being a periodic 
function p(t); this is the case for certain values of a con- 
stant external forcing / p7|. or when the neuron is cou- 
pled to another periodic spiking or bursting neuron. So 
in a next step we consider the fast (w, u!)-subsystem with 
periodic forcing 

v(t) = w - v 3 + 3« 2 + p (t) 
w(t) = 1.8-5v 2 -w, (10) 

or, in Lienard form 

p(t)=v + f(v)v + g(v), (11) 

where / and g are defined as above. According to |29| . 
theorem 4.3.1 and 4.3.3, equation Ijlll) has a non-constant 
periodic solution with the same period as the forcing term 
p{t) if several conditions (i-vi) (given in Appendix IB 4|) 
are fulfilled. Hence, we make the following proposition: 
When p(t) is a T-periodic function, then l(TTJ> has a peri- 
odic solution with the same period T. The detailed proof 
is in Appendix IB 41 

For illustration, we show the behavior of the 
Hindmarsh-Rose equations JBHSJ) excited with a constant 
forcing I(t) — I (see Fig. 0), thus again z(t) is a peri- 
odic function. This is an example for a typical intrin- 
sic burs ter, where the slow variable slaves the fast sub- 
system [23. The theorem presented above proves that 
the fast (v, w)-subsystem of the Hindmarsh-Rose model 
shows mode locking when it is driven by a periodic forc- 
ing, no matter what frequency or amplitude this forcing 
has. 

If we assume e to be small (~ 0.006 as used by 32]), 
z can be treated like a slowly varying parameter in the 
(v,w) system ()10|l . To understand the influence of v on 
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FIG. 7: (Color online) v(t) (black) and z{t) (red/gray) with 
I = 1.0. When z(t) is periodic, v(t) oscillates with the same 
period; thus the slow variable z(t) slaves v(t). 
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FIG. 8: (Color online) v(t), z(t) and I(t). The system shows 
a post-inhibitory rebound burst after a hyperpolarizing step 
of I(t) = —0.5 and a duration of 70ms. 
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we write it in the form of a relaxation pro- 



-e(z - 



>(«)), 



(12) 



where Zoa(v) = 3.3(t> + 1.56). As e is very small, the fast 
dynamics of u has only little influence on equation l|12Jl . 
as its own dynamics is too slow to follow the fast z 00(11). 
Because of this, z will arrive in a steady state if v is 
spiking quickly, and then just oscillate with a very little 
amplitude around it. This effect leads to the so-called 
spike frequency adaption (for details see Q), which is 
also observed in biological neurons [24| . For small enough 
values of I(t) this effect can also lead to an effect similar 
to intrinsic bursting [27l |. 



Post-inhibitory rebound bursts in the 
Hindmarsh-Rose neuron 



low threshold calcium current It (see Fig. |2Jj) Q , so the 
variable z can be interpreted as a calcium current. 




inhibitory 
current: 



Tc cell 



Re cell 



input: 




Post-inhibitory rebound bursts are a dynamical fea- 
ture being also present in the biophysical neuron model 
Q • As it is crucial for spindle oscillations and for the in- 
formation transfer, we will describe it in detail here, see 
also Fig. [S] If the Hindmarsh-Rose model 0] is hyper- 
polarized for a period similar to the burst duration, the 
adaption current z will decrease below its steady state. 
If the hyperpolarizing current is released, the z current 
stays below its steady state for some time due to its slow 
timescale. As a consequence, the model behaves as if 
some extra current had been applied, so it will cross the 
saddle node bifurcation in Fig. H3 and stay on the limit 
cycle. Due to the spikes, z will increase above its steady 
state, stopping the burst (see Figs. l8l9[) . Due to the slow 
timescale, a postinhibitory rebound burst will only ap- 
pear if slowly-decaying dynamic synapses are used, as 
otherwise z will not decrease below its steady state. In 
the biophysical model, the same effect occurs due to the 



FIG. 9: (Color online) Right: Reciprocally coupled RE and 
TC cell. Left: If the excitor cell gets activated by an ex- 
ternal pulse, it excites the inhibitor cell, which leads to an 
inhibitory current. This current hyperpolarizes the excitor 
cell and evokes a rebound burst in it, and the mechanism 
repeats, which leads to self- sustained oscillations. 



In reciprocally coupled neurons, the effect described 
above leads to self-sustained oscillations if the excitor 
gets activated by a single input spike; this effect is an 
example for so-called hard-excited self-oscillations |33| . 
During this oscillation, the two neurons exhibit anti- 
phase synchronization. We show this for the case of two 
Hindmarsh-Rose neurons, which are coupled like in a typ- 
ical thalamocortical circuit (see Fig. [§J| . 
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FIG. 10: (Color online) Reciprocally coupled inhibitor and 
excitor; without the h current, the oscillation does not termi- 
nate. 



D. Motivation of a simplified calcium current: 
Extending the Hindmarsh-Rose model 



The self-sustained oscillations in Fig. ^| are stable and 
do not terminate; this is contrary to the experimental ob- 
servations of Le Masson et al. @ . A spindle oscillation or 
a burst consists of a series of metastable self-oscillations 
enhanced by rebound bursts which terminate after a few 
seconds. Above we argumented that the slow variable z 
can be interpreted as a calcium current, and such cor- 
responds to the low threshold calcium current It in the 
biophysical model. The schematic diagram Fig. [S] shows 
that the reason for the termination of the oscillation is 
a further ion current K H II M- As the original 
Hindmarsh-Rose model does not incorporate a variable 
which corresponds to the Ih current in the biophysical 
model, we will extend the Hindmarsh-Rose model by an 
equation motivated directly by the dynamics of Ih in the 
biophysical model. Destexhe et al. [4( model Ih by a 
double activation kinetic, consisting of slow and fast ac- 



Fa - 



* Si 

Fx, 



(13) 
(14) 



were So/i and f /i represent the closed and open states 
of the slow and fast activation gates, respectively. The 
open state gates are assumed to have n bindings for Ca 2+ 
which lead to the open bounded gates £2 and F 2 , 



Si 
Fi- 



nCa 2+ 
nCa 2+ *- 



~>S 2 
F 2 . 



(15) 
(16) 



In this model the activation function of Ih is shifted dur- 
ing the oscillatory phase by the entry of Ca 2+ , and thus 
terminates the oscillation. Destexhe et al. £j find that 
the length of the silent phase and of the oscillatory phase 
were directly proportional to the time constant of intra- 
cellular Ca 2+ binding to Ih channels, k 2 . Further it is 
assumed that the binding of Ca 2+ is critical for the on- 
set and termination of the oscillatory phase So the 
length of the oscillatory phase depends on the rate of rise 
of the variables S 2 and F 2 , whilst the length of the silent 
phase depends on the rate of relaxation of S 2 and F 2 back 
to their resting values (see Fig. 1111 and Q). According 
to (IA20I) . both the length of rising and relaxation of S 2 



(red) and F 2 are proportional to fcT • As Fig. 1111 shows. 
F 2 only displays small variations of amplitude compared 
to S 2 and therefore plays a less important rule. In our 
further simplification we will neglect the influence of F 2 . 
The slow variable S 2 slaves the system and switches be- 
tween the oscillatory and the resting state leading to wax- 
ing and wanning oscillations 4]. Our simplified model 
will be motivated from the kinetic equation of S 2 . As 
the transition from So to Si is much faster than the one 
from S\ to S 2 , we will assume it to be instantaneous, 
what leads to S 2 ~ Sq, so we get 



So + nC a 



2+ 



S 2 



(17) 



In equation (5) in |j] , the number of binding sites n for 
Ca 2+ is assumed to be 2, in our model we use n = 1 with 
the backward rate k 2 = 4 * 10~ 4 , the forward rate k\ = 



and C 



[Ca, 



- 10 ^r aim >_y = g^ 10 _4 , where [Co m ] is the intracellular 
calcium concentration, thus altogether we get 

S 2 = -k 2 [S 2 - CS ] . (18) 

The essential features of equation (|T%)l can be summa- 
rized as follows: S 2 opens and closes proportional to the 
same rate constant k 2 ; due to C, the activation depends 
critically on the concentration of intracellular Ca 2+ . This 
features should also be present in our simplified model. In 
the Hindmarsh-Rose model the calcium current is mim- 
icked by the variable z [l4| , in our model for the h current 
will get activated by the variable z. The rate of activa- 
tion and deactivation will be equal to the constant used 
in the biophysical model. For convenience and to be con- 
sistent with the notation used in biophysics, we call this 
current h. 



t-tc 



-k 2 {h 



TC 



0. 



9- z T c))- 



(19) 



8 




6000 



8000 



10000 

t (ms) 



12000 14000 



FIG. 11: (Color online) Top: Time course ol the gating vari- 
ables S2 (red) and F2 (green) during a spindle cycle. Bottom: 
The membrane potential Vtc of the biophysical model 
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FIG. 12: (Color online) Top: Time course of h during a spin- 
dle cycle. Bottom: The Membrane potential vt c of the ex- 
tended Hindmarsh-Rose model. 



As z decreases when the cell gets hyperpolarized, we use 
the difference between the maximum magnitude of z and 
z itself as an activation term in 119|) . In order to repro- 
duce the behavior of the biological TC cell, we extend 
the Hindmarsh-Rose equations as follows. 



vtc 
wtc 
ztc 

h-TC 



3vi 



wtc — v TC t o v tc 
1.8 — 5v TC — wtc 
e [3.3(«rc + 1-56) - 
-0.0004(/i TC + 0.8; 



ZTC 



l>TC 



ztc] 

1(0.9 -z)). 



(20) 
(21) 
(22) 
(23) 



Equations i|2QI23[l define our central model, and will be 
used to model the reciprocally coupled neurons in the 
remainder. As Fig. ll2l shows. the time course of the vari- 
able h is very similar to the time course of S2 in the 
biophysical model. To understand the influence of the 
variable h we proceed in the same way as before. The 
timescale of the dynamics of h is less than a tenth than 
the timescale of z. As declared above, the reason for the 
self-sustained oscillations is the after-hyperpolarization 
activation by the variable z. In the presence of h, the 
slowly varying parameter in the (v, w) system is the sum 



(z(t) + h(t)). As h(t) gets activated by z(t), after some 
rebound bursts the sum (z(t) + h(t)) is below the thresh- 
old of the (v, w) system, as a consequence, the oscillation 
terminates. When the excitor cell is inactive, there are 
no more inhibitory currents, that means that t is inactive 
too, so h decays slowly until it is small enough that the 
the system can get activated again by an input spike. 



E. Reciprocally coupled TC-RE neurons using the 
extended Hindmarsh-Rose model 



As the h current is absent in the biological RE neuron, 
we have to extend the Hindmarsh-Rose system by the h 
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equation ill'MI) for the TC neuron only. 



Vtc 
wtc 

ZTC 

hrc 
vre 

Wre 

zre 



wtc 



3v 



TC 



ztc — h-Tc + Igaba + I(t) 



1.8 — 5v TC — wtc 
0.006(4(>tc + 1-56) - ztc) 
-0.0004(/i TC + 0.88(0.9 - z)) 



WRE 



J RE 



3v 



RE 



zre + Iglu 



1.8 - 5v RE - wre 
0.006(4(>fl E + 1.56) - z RE ) 



(24) 



Igaba and Iglu are governed by equations l|B4(l and 
(|B5|) in appendix IB II respectively. The external in- 
put was modeled as described above for the biophys- 
ical model, with modified parameters adapted to the 
Hindmarsh-Rose model, f was 1/100 and To was 30ms. If 
we stimulate our model by such a modified Poisson pro- 
cess, we have a computational model for the experiment 
of Le Masson et al. 2} , which allows us to compare our 
computational model with the experiment and the bio- 
physical model Q , so that we can validate our extension 
of the Hindmarsh-Rose model. 

A comparison of Fig.HSand (@, Fig. 2) shows that our 
system reproduces the experimental results of Le Masson 
et al. quite well. During the spindling state, the firing 
pattern of the TC cells (excitor), which is very different 
from the input, shows that the information transfer of 
the input is low. The question whether this low transfer 
is still reliable is calculated in the same way as for the 
biophysical system. With a strong inhibitory feedback, 
the Tsn was low, showing that most of the TC spikes 
where not triggered by an input spike and thus that the 
system is not transferring spikes in a one-to-one manner. 
This result encouraged us to screen the strength of inhi- 
bition to test if recurrent feedback inhibition could be a 
way to control the precision of input-output-transfer in 
a wide range, or if there is just a switching between a 
relay state and masquerading state. Our numerical in- 
vestigation leads to the following results. An increase of 
the inhibitory coupling strength leads to a smooth de- 
crease in Tsn from a maximum value of around 1 to a 
minimum value of less than 0.3 (Fig. I14JI . Despite this 
significant decrease in the reliability of spike transfer, the 
efficiency of spike transfer Tte was not significantly di- 
minished by the strength of the inhibitory coupling (Fig. 
I14JI . Thus inhibitory feedback has a direct decorrelating 
effect, which is able to reduce the reliability of the spike 
transfer. 



Dynamical behavior of coupled TC-RE neurons 
using the extended Hindmarsh-Rose model 



As we argued in llll Al for the Hindmarsh-Rose model, 
Destexhe et al. Q use singular approximation to char- 
acterize the spindle oscillations as a transition between 
a hyperpolarizing stationary phase and an oscillatory 
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FIG. 13: (Color online) Spontaneous spindle activity in our 
computational model (the numerical values of the voltages are 
scaled by a factor of 30 to match with the biophysical values) . 
Similar rescalings have been also necessary in the underly- 
ing Hindmarsh-Rose model (see, e.g. |2^|). a) The computa- 
tional circuit: A TC cell (excitor) modeled by the extended 
Hindmarsh-Rose model is reciprocally coupled to a model RE 
cell (inhibitor) represented by the original Hindmarsh-Rose 
model. The TC cell receives artificial synaptic retinal bom- 
bardment modeled by a Poisson distributed spike train. Like 
in the experiment the system shows spindle activity, b) Detail 
of a), c) Like in Fig. GO 



phase. Here we will apply this method to our simpli- 
fied model of a thalamocortical oscillator. In Fig. the 
h current evolves according to a much slower time scale 
as the dynamics of the coupled TC-RE neurons. 

We assume h to be a slow varying parameter in 124|) 
with I(t) = 0, so we consider the limit — * in 11911 . 
We will consider h as bifurcation parameter of l|24|) in an 
interval between -1 and 1.5 as a parameter. In this range 
of h the dynamical state of the system can be divided 
in three areas: For high values of h the system is in a 
stable resting state, for intermediate values of h both a 
stable resting state and a stable oscillatory state exist, 
separated by an unstable limit cycle. For negative val- 
ues of h the stable resting state gets unstable, while the 
stable oscillatory state persists. The transition between 
the stable fixed point and the oscillatory state occurs by 
a subcritical Hopf bifurcation [13j . The existence of a re- 
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FIG. 14: (Color online) Percentage of output spikes triggered 
by an input spike. An increase of the strength of the in- 
hibitory coupling gcABA comes along with a significant de- 
crease in the reliability of the information transfer. As in 
the experiment of Q and the in the biophysical model, the 
transfer efficiency does not vary significantly. 



gion of parameters where stable solutions overlap is typ- 
ical for a subcritical Hopf bifurcations; our system shows 
this typical behavior in a wide range of the bifurcation 
parameter. In this states the system is bistable, i.e. a sta- 
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FIG. 15: Bifurcation diagram of Il21t with h treated as 
a slowly varying parameter; solid lines indicate stable fixed 
points, dashed lines unstable fixed points, filled circles stable 
limit cycles, open circles unstable limit cycles. For extreme 
values of h, the system is monostable, i.e. either a stable fixed 
point or a stable limit cycle coexist. 

ble fixed point and a stable limit-cycle coexist (see Fig. 
I15p . In this area, the state of the system depends on its 
history: If the initial conditions or the past state of the 
system lie within the basin of attraction of the limit cy- 
cle, the system exhibits stable oscillations; if the starting 



point lies within the basin of attraction of the fixed point, 
the system rests on the fixed point. These two attractive 
areas are separated by an unstable limit cycle. Further, 
the system can be switched between these two states by 
an input pulse. 



G. Discussion of the transfer-properties of the 
extended Hindmarsh-Rose model 

The h current switches the system between two exist- 
ing dynamical states, namely a fixed point and a bistable 
state where a limit cycle and a fixed point coexist. With- 
out this further current (i.e. h = 0), two reciprocally 
coupled Hindmarsh-Rose neurons exhibit so-called hard 
self-excitation [3^ as response to an input spike, result- 
ing in stable self-sustained oscillations. The h current 
makes this oscillatory state metastable, resulting in wax- 
ing and waning oscillations. An increase of the inhibition 
changes the responsivity of the system: Without, or with 
only little inhibition, the system responds to an input 
pulse with an output spike, so that spikes get transmit- 
ted in a one-to-one manner. With high gain in inhibitory 
feedback, the system responds to an input spike by an 
output burst (range B in Fig ll5|l . followed by a silent pe- 
riod where transmission is totally blocked (range A in Fig 
I15|l . During the burst, the system is in an autonomous 
self-oscillatory state which masquerades the input, this 
leads to an increase of the T$n- As the system is still 
excitable, but with a different responsivity the Tte does 
not change as strong as the T$n- So the transmission 
behavior depends critically on the rate of rise and fall of 
the h equation. This insight might also help to improve 
the biophysical model in terms of information transfer. 



IV. CONCLUSIONS 

We have shown that an extended Hindmarsh-Rose 
model is able to reproduce the behavior of a biological 
thalamocortical relay neuron in the recent experiment 
of Le Masson et al. To gain more insight into the dy- 
namical mechanisms, a simplification or reduction of the 
detailed biophysical model, which did not convincingly 
reproduce the experimentally observed decrease of the 
signal-to-noise ratio, was necessary. The widely used 
Hindmarsh-Rose model however does not show the char- 
acteristic waxing and waning oscillations. Further it does 
not exhibit the quiescence periods necessary for transfer 
of information. Especially, we have analytically proven 
that the fast subsystem of the Hindmarsh-Rose model 
gets slaved by periodic forcing, which can lead to intrin- 
sic bursting if the forcing is slow. To account for the 
waxing and waning mechanism, we propose an extended 
Hindmarsh-Rose neuron model for the TC cell, directly 
motivated from the biophysical model, taking the low 
threshold calcium current explicitcly into account. 

In this paper, we have demonstrated that our extended 
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Hindmarsh-Rose model serves as a computational model 
for the setup in the Le Masson et al. experiment, show- 
ing that the information transfer can be adjusted within 
a wide range by the gain of the recurrent feedback in- 
hibition. Our numerical investigations strongly suggest 
that the low threshold calcium current plays an impor- 
tant role for the information transfer in thalamocortical 
circuits. From the technical simplicity of our approach 
and its agreement with the experiment, this approach 
may give rise to similar models for neural systems where 
a third timescale is apparent, and a second slow degree 
of freedom is necessary to describe the dynamics. 
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APPENDIX A: THALAMOCORTICAL CIRCUIT 

The thalamocortical circuit consists of a pair of TC and 
thalamic RE neurons connected as shown in Fig^ For 
each TC and RE cell several ion currents were included 
which will be described in detail here. All intrinsic cur- 
rents are described by the same general equation [20l |: 



g j m M h N (V-E j ). 



(Al) 



Here the current Ij is the product of the maximal con- 
ductance, gj, activation m and inactivation h variable, 
and the difference between the membrane potential V 
and the Nernst reversal potential Ej . The powers N and 
M are the respective numbers of ion channels to be open 
synchroneously. The gating of a membrane channel is 
described by a first order kinetic scheme 



c o 

c m. o. 



(A2) 
(A3) 



where O and C stand for the open and the closed state of 
the gate, with a(V) and j3(V) as the transition rates. Ac- 
cording to the Hodgkin-Huxley model |2jj , the variables 
m and h represent the fraction of independent gates in 
the open state and are described by simple first order 
differential equations, 



1 



T m (V) 
1 



(m - m x (y)) 
(h-hooiV)). 



(A4) 
(A5) 



The steady states m^V) and hoo(V) and the time con- 
stants T m /h(V), respectively, can be written as functions 
of the transition rates a and f3, using x S {m, h}, as 



Tx{Y) a x (V)+p x (V) 
Xoo(V) = a x {V)r x {V). 



(A6) 



The TC cell 



Here we summarize the results of the previously men- 
tioned investigations of TC cell membrane properties. 
The membrane potential of the TC neuron is given by 

C m Vr = —It l — It — hi — Ir Na — It k 

^GABA„ - -^GABAb • ( A 7) 

The area of a TC cell membrane is about 2.9 * 10 -4 cm 2 
which is according to a cell capacity of C m — 0.2%nF. 
All constants used in the simulations correspond to a cell 
of this size. 

The leakage current It l is the only passive current, 
and is governed by the Ohm law 



It l = 5t l (^t - E L ) 



(A8) 



where gx L = Q.05mS', El = —86mV. 

All the other currents are active currents, and more 
complicated, as detailed below. 

The sodium current I^ a has the form 



g Na m{tfh{t){V T - E Na ), 



(A9) 



with gT Na — 30mS and E^ a — 50mV. In addition we 
have for the gating variables m(t) and h(t) 



hoo(V T ) = 
r m (V T ) - 
T h (V T ) = 



Cim{V T 



where 



a m (V T ) = 0.32 



a m (V T ) + m (V T ) 

a h (V T ) + MVt) 
1 

a m (V T ) + Pm(V T ) 
1 

a h {V T ) + h (V T y 
(V T + 37) 



(A10) 



l-exp(-™) 



(All) 



(V T + 10) 



a h {V T ) = 0.128 exp(- Vt + 33 ) 

5 

h (V T ) 



exp(- 



1 



With equations l|A4IA5|l the sodium current is completely 
described. 
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The potassium current Ik has the form 



It k =g TK m(t) 4 (V T -E k 



(A12) 



with gx K — 2.mS and E k — — 95mV. The gating and 
transition variables are given by 



moo 
r m (V T ) 
ot m (V T ) 



oi m {V T ) 



a m {V T ) + P m (V T ) 
1 

a m {V T ) + f3 m (V T ) 
(V T + 35) 



(A13) 



0.032 



l-exp(-i^5) 



0.5exp( 



V T + 40, 



With equation (| A4|) the potassium current is completely 
described. 

The low threshold Ca 2+ current It is taken to be 



It = gc a m(tfh(t){V - E t ), 



(A14) 



where gc a — 1.75mS and the reversal potential Et de- 
pends on the Ca 2+ concentration inside ([Ca]i n ) and out- 
side ([Ca] out ) the cell. It is defined by the Nernst equa- 
tion 



E T = 1000*:^ In ^2*. (A15) 



where R = 8.31441 



K-mol ' 



T = 309.15if,F = 96489^, 
and [Ca] ouf = 2mM is considered to be constant. In the 
remainder, all concentrations will be denoted by square 
brackets, following an usual convention (see e.g. [l5|l- 
The calcium dynamics in the cell are described by a sim- 
ple model which was introduced in ^3] , 



In addition, for this current the constants are 
,Fr + 65 x 



TOoo(Vt) 
T m (V T ) 

ai(V T ) 

K{V T ) 
a 2 (V T ) 
t 2 (V t ) 



1 +exp(- 



7.. 



0.15 moo (V T ) ( 1.7 + exp(- VT +f ~ 8 ) 



13. 



exp 



{V T + 162.3) 
17^8 



/0.26 



'0.25 + exp 
1 



(V T + 85.5) 
6.3 



0.5 



r 2 {V T ) (K(V t ) + 1) 
62.4 



1 



exp(^) 



(A18) 



The numerical values are fits to experimental data |18j . 

The hyperpolarization- activated cation current Ih fi- 
nally is described by 



h = 9h(Si + 5 2 )(Fi + F 2 )(V T - E h 



(A19) 



where gt — Q.15mS and Eh = — 43mV. The gating 
variables for this current are governed by 

Si = a s {Vt)S -/3 s (Vt)S 1 + k 2 (S 2 -CS 1 ) 
F 1 = a f (V t )F -f3f(V t )F 1 + k 2 (F 2 -CF 1 ) 

s 2 = -ms 2 -cSi) 

F 2 = -k 2 {F 2 -CF x ), (A20) 



where C 



[Ca,. 



4 * 10 4 , the rate constants 



5*10" 4 ' ■ 

a s /f(Vt) and (i s /f(Vt) are related to the activation func- 
tion and the time constants t s / f 



a s (V t ) = Hoo/ts 
a f (V t ) = H^/tj 
P s (V t ) = (l-JZooVr, 



where 



d[Ca] 
df 



-Ah 



K T [Ca] in /([Ca] in + K d ), (A16) H ~ = + cxp ( 6 . 5 J J 



(A21) 



(A22) 



where A = 0.179 



mmol 



K 7 



10 



-4 mmol 



and Kri = Kq 



t s = cxp((y t + 183.6)/15.24) 
'V t + 158.6 



ms-fiA 5 

The gating variables are governed by the following system 
which was proposed by Wang |^| 



rh = (r m (Vr)) _1 (m - TOoo(Vt)) 
/i = a 1 (V T )(l-/i-d-/s:(V T )/i) 
d = q;2(Vt) (-fr(F T )(l- h-d) - d) . (A17) 



exp 



11.2 



1 + exp 



Vt+75 
5.5 



This makes the biophysical description of TC cell 
complete. While it contains a lot of simplifications, 
for example the single compartment assumption, it 
reproduces almost all important features of TC cells 
quite good. We will use this model in our computer ex- 
periment to verify theories obtained by using simplified 
models. For more details see [l7|. 
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2. The RE cell 

The membrane potential of the thalamic RE-ncuron is 
governed by the cable equation: 

C m Vii = —Ir l — Its — Ik[Co\ — Ican — Ir No . — Ir k 
-Iglu- (A23) 

For the RE cell the membrane capacity is C m = 0.143nF, 
the leakage conductance is g Rl — 0.05mS and the reversal 
potential of the leakage current is El = — 80mV. This 
model of thalamic RE cell was introduced by Destexhe 
et al. in [Sfij . 

The equations for the sodium Ir Nu and potassium Ir k 
current in the RE cell are the same as for the TC cell, 
except that the conductance gNa — lOOmS* is differ- 
ent. The description of the RE cell potassium current 
is reached by putting gx — lOmS. 

The kinetics of the low threshold Ca 2+ current It s 
was established as model for voltage clamp data on rat 
RE cells using a so-called m 2 h formalism 



Its = g Ts m(t) 2 h(t)(V R - E Ts ), 



(A24) 



where gx s — l-75fiS and Et s depends on the calcium 
concentration in the same way as described for the TC 
cell. The calcium dynamics are described by equation 
(|A9(I . when It is replaced by Its- In addition we have 
for this current 



T-oo{V R ) = ( 
T m {V R ) = 1 



1 + exp( 
1 



exp( 



hUV R ) = l+exp( 



Vr + 52 . 

7.4 
Vr + 27 . 
10 ' 
Vr + 80^- 1 



(A25) 



exp(- 



Vr + 102 
15 



85 1/ Vr + 407 \-i 

The full kinetics of the It current is reached if m,oo, T m , 
/loo and Th are inserted in equations (|A4IA5|1 . 

The calcium dependent currents: The RE cell pos- 
sesses further two calcium dependent currents the slow 
K + current Ik[Co,]i and the slow nonspecific cation cur- 
rent Ican- According to fl?) they are modeled as 
voltage-independent currents described by equation (|A1() 
with M=2 and N=0. The activation variable m obeys 



m = a[Ccii] 2 (l — m) — (3m, 



(A26) 
the 



where a and (3 are rate constants and [Co^j is 
intracellular calcium concentration. For Ik[Co\ 
48m,s' 1 m,M~ 2 and f3 = 0.03ms" 1 , for I C an 

With this the 



a 



a = 20ms~ 1 mM- 2 and (3 = 0.002ms" 1 
description of the RE cell is complete. 

3. Synaptic currents 



Kinetic models of synaptic currents have been di- 
rectly fitted to experimental data 0, based on whole 



cell recorded synaptic currents obtained in hippocampal 
neurons [33 • The GABA A and GLU (=glutamate) re- 
ceptor mediated currents are both represented by a first 
order kinetic scheme 1381 



C + T^O 



(A27) 



where the transition between closed (C) and open (O) 
states depend on the binding of the transmitter (T). The 
current is given by 



g S yn[0}(t)(V(t)-E syn ), 



(A28) 



where g syn is the maximal conductivity and E syn is the 
reversal potential. For AMPA receptors E syn = OmV 
and for GABAa receptors E syn = — 80mV. g syn = 0.05 
for GABAa and g syn = 0.1 for glutamate synapses. 
[O] (t) is the fraction of open channels 



d[Q](t) 
dt 



= a{l-[0}(t)}[T](t)-P[0](t) ) (A29) 



and [T] (t) is the concentration of transmitter released 
from time t to time t max , 



[T]{t) = AQ(t r , 



t)Q(t), 



(A30) 



where 0(x) is the Heaviside function. The synaptic 
parameter values used in the simulation are A = 0.5 
and t max = 3ms for GABAa and 9ms for AMPA 
synapses, the rate constants were chosen as a — 5ms and 
(3 = 0.16ms for GABAa synapses and a = 0.94ms and 
(3 = 0.18ms for GLU synapses. The activation scheme 
of GABAg receptors is more complex, as it involves the 
activation of K + channels by G proteins (for details see 
|39jl- The model used by Le Masson et al. is a mod- 
ified version of the GAB A b kinetic model introduced in 

m 



Ro+T 
R + Go 
G 

Ci +4G 



R^D 

RG — ► R + G 

G 

O, 



(A31) 
(A32) 
(A33) 
(A34) 



where the transmitter, T, binds to the receptor, Rq, lead- 
ing to its activated form, R, and desensitized form, D. 
The G protein is transformed from an inactive (GDP- 
bound) form, Go to an activated form, G, catalyzed by 
R. Finally, G binds to open the K + channel, with four 
independent binding sites. With some assumptions (see 
0) the kinetic model for this system reduces to 



IgabAb 

d[R] 
~dT 
d[G] 
dt 



gGABA B - 



[GY 



[G] 4 + K d 
Ki{l-[R])[T\-K 2 [R] 

K 3 [R]-K 4 [G], 



[V(t)-E k ] (A35) 



where [R] is the fraction of activated receptors and [G] 
is the concentration of G proteins. In these equations 
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9gaba b , Ki = 0.52, K 2 = 0.0013, K 3 , K 4 = 0.033 and 
K d = 100. 

The strengths of the synaptic couplings gGABA A/B and 
9am pa are varied in different simulations. 

APPENDIX B: THE REDUCED MODEL 



the reversal potential of the GABAa synapse was 
rescaled by a factor of 30 to adapt it to the scale of the 
Hindmarsh-Rose model. 



2. Transformation to a Lienard system 



1. Synaptic Currents: 



We begin with the autonomous (v,w) subsystem (|6I7[) . 



The simplest way to model synaptic currents is a two 
state scheme of the binding of a neurotransmitter T to 
postsynaptic receptors |3| . So the receptor mediated cur- 
rents are both represented by a first order kinetic scheme 

C + T^O, 

where the transition between closed (C) and open (O) 
states depend on the binding of the transmitter (T) with 
forward and backward rates a and (3 respectively. As in 
[30| we only consider voltage- independent rate constants. 
With this assumptions we get the following kinetic equa- 
tion 

d\Q] 
dt 



a{l-[0])Tr{t)-(3[G] = aTr{t)-[0]aTr(t)-(3[0}. 

(Bl) 



The synaptic current is given by 



g syn [0](t)(V(t)-E syn), 



(B2) 



where g syn is the maximal conductivity and E syn is the 



reversal potential. For glutamate receptors E s 



OmV 



and for GABAa receptors E s 



-50mV. The release 



and clearance of transmitter are extremely fast processes 
compared to the open/close kinetics, resulting in a very 
brief presence of transmitter in the synaptic cleft |31j . 
which allows us to consider O to be constant during the 
transmitter release. For simplicity we assume the trans- 
mitter release time course to be a square pulse AQ(V pr e) 
which occurs when the presynaptic neuron fires a spike, 
i.e. when the presynaptic potential V pre gets depolarized 
towards positive potentials. If we assume that initially 
all synaptic channels are in the closed state, we get 



d[0] 
dt 



jaAQiVpre) - P[0], 



(B3) 



were 7 = a A. As GAB As synapses contribute less than 
a tenth of the total inhibitory synaptic current, they 
are neglected here. The rate constants and amplitude 
were taken from the biophysical model IA 31 The con- 
stants for the GAB A synapse were modified because of 
the absence of the GABAb synapse: a — 5ms -1 and 
(3 = 0.05ms -1 . For the GLU synapse we used the same 
constants as in the biophysical model: a = 0.94ms -1 
and /3 = 0.18ms -1 . We chose A = 0.5, and respectively 
7 = 2.5 for GAB A synapses and 7 = 0.47 for glutamate 
synapses. The synaptic currents are governed by the fol- 
lowing equations: 

Igaba = 9gaba[0]gaba{vtc - 2.5) (B4) 
I glu = 9glu[0]glijvre, (B5) 



v(t) = w-v 3 + 3v 2 
w[t) = 1.8 -5v 2 - w. 



(B6) 
(B7) 



Then v — w + (6t> — 3v 2 )v, w — v — (6i> — 3v 2 )v, and 

w = -v + (6v - 3v 2 )i) + 1.8 - bv. (B8) 
Inserting l|B8|l in (|B6|I yields 



= i + (1.8-6v + 3v 2 )v + (v 3 + 2v 2 -1.8) 
^0 = v + f(v)v + g(v). (B9) 



(IB9|) is the so called Lienard form of equations l|B6(l and 
(fB7|) . where f(v) = l-6v + 3v 2 and g(v) = u 3 +2u 2 -1.8. 

3. Stability of the equilibrium points 

The stability of equilibrium points in the fast (v, w)- 
subsystem is investigated by a linear approximation to 
the system (|6l7ll . Suppose the equilibrium point has the 
w-coordinatc vq, then the linear approximation is: 



(B10) 



where (x, y) are new coordinates whose origin is in the 
equilibrium point, and A(vq) is the Jacobian in vq 



(BID 



The kind of equilibrium point may be determined by the 
signs of the trace Tr(A(uo)) = — 3v 2 + 6vq — 1 and deter- 
minant Det(j4(w )) = 3v 2 + 4v Olil|. The eigenvalues 
A1/2 of the Jacobian A(vo) are expressed by 

A1/2 = l{Tr(A(v a )) ± v/(Tr(A(«o))) 2 -4Det(A( Uo ))} 

(B12) 

Tabled gives the type of equilibrium point according 
to the region vo belongs to; here Tr_ = (3 — \/6)/3 and 
Tr + = (3 + V6) /3 are the negative and postive zeroes of 
Ti(A(vq)), respectively. 



4. Periodic solutions for periodic forcing 

The conditions for constants b, m, M > are 
(i) for |u| > b,f(v) > m 
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TABLE I: Equilibrium points depending on vo (see text). 
Only stable nodes and unstable spirals and saddle nodes oc- 
cur. 



region 


values of vo 


sign of 


sign of 


type of 






TvA(v ) 


DetA(vo) 


node 


I 


vo < -4/3 




+ 


stable 


II 


-4/3 < v < Tr_ 






saddle 


III 


Tr_ < vo < 


+ 




saddle 


IV 


< v < Tr+ 


+ 


+ 


unstable spiral 


V 


Tr+ < vo 


+ 




saddle 



(ii) V»e R,f(v) > -M 

(iii) for \v\ > b,vg(v) > 

(iv) g(v) is monotone increasing in (— oo, —6) 

(v) \g(v)\ — ► oo for \v\ — > oo 



(vi) g{v)/G{v) — ► for |w| — > oo, 
where G(w) = J g(u)du. 

The proof can be sketched as follows: 

(a) The existence of M is easy to establish as f(v) has 
a minimum at v = 1, further as /(u) is a convex 
function this minimum is the global one. If a 
M > |/(1) = -2| is chosen (ii) holds. 

(b) Set M = m, then 6i may be chosen arbitrarily. 

(c) Next vg(v) is a fourth degree polynomial where the 
coefficient of the w 4 -term is positive. So there exists 
a 6 2 with vg(v) > if v > b 2 > 0. 

(d) As g(v) is a third degree polynomial, there is a 
63, such that g(v) is monotonously increasing in 
(-00, -63) and (63, 00). 

Choose b = max(6i, b 2 , 63) then (i), (iii) and (iv) hold. 
Finally (v) and (vi) are obvious; so the proof is complete. 
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